Spatial organization in cyclic Lotka-Volterra systems 
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We study the evolution of a system of AT interacting species which mimics the dynamics of a 
cyclic food chain. On a one-dimensional lattice with N < 5 species, spatial inhomogeneities develop 
spontaneously in initially homogeneous systems. The arising spatial patterns form a mosaic of 
single-species domains with algebraically growing size, £(t) ~ t a , where a = 3/4 (1/2) and 1/3 
for N = 3 with sequential (parallel) dynamics and N = 4, respectively. The domain distribution 
also exhibits a self-similar spatial structure which is characterized by an additional length scale, 
C(t) ~ t 13 , with (3 = 1 and 2/3 for N = 3 and 4, respectively. For TV > 5, the system quickly reaches 
a frozen state with non interacting neighboring species. We investigate the time distribution of the 
number of mutations of a site using scaling arguments as well as an exact solution for TV = 3. Some 
possible extensions of the system are analyzed. 



PACS numbers: 02.50.Ga, 05.70.Ln, 05.40.+j 
I. INTRODUCTION 

The classic Lotka-Volterra equations mimic the 

dynamics of interacting species such as predator-prey 
systems. These equations are rather successful in pre- 
dicting density oscillations which are known to exist in 
Nature. For spatially inhomogeneous situations, Lotka- 
Volterra equations |^| are straightforwardly generalized 
to diffusion- reaction equations |3|; these equations were 
widely applied to more complex ecological processes. 
However, such an approach ignores spatial correlations 
and therefore fails to predict the development of spa- 
tial heterogeneities in initially homogeneous systems. For 
chemical processes, the crucial role that spatial hetero- 
geneities play in governing the kinetics has been appre- 
ciated over the past decade, see e.g. Q and references 
therein. Therefore, in low spatial dimensions the mean- 
field like rate equations approach (analog of the Lotka- 
Volterra equations in chemical kinetics) fails to provide 
the correct asymptotic behavior. Indeed, a homogeneous 
initial state evolves to a strongly heterogeneous state, 
namely to a coarsening mosaic of reactants which con- 
fines the actual microscopic reaction to the interfacial 
regions between domains, and therefore the kinetics are 
significantly slowed down. Similar spatial organization 
was recently reported in Lotka-Volterra systems pT| . 
However, theoretical understanding of these systems is 
still incomplete. 

In this study, we consider the evolution of an N species 
food chain, where every species plays the role of prey and 
predator simultaneously. The food chain is thus assumed 
to be cyclic; e.g., in the 3-species system, A eats B, B 
eats C, and C eats A. Every "eating" event leads to 
duplication of the winner and elimination of the loser, 
therefore the 3-species food chain is symbolized by the 
reaction scheme 

A + B^2A, B + C^2B, C + A -> 2C. (1) 

The corresponding stochastic process is well defined on 



a lattice, where the interaction is restricted to nearest 
neighbor sites. Initially, every lattice site is assumed to 
be occupied; clearly, the lattice then remains fully occu- 
pied. 

Given the simplicity of the reaction process (|l|), one 
anticipates that it can provide a caricature description 
of a number of phenomena in Nature and Society. One 
example is the voter model |lj,[l3|], which is applica- 
ble to chemical reactions on catalytic surfaces [l3,ff5l. 
This model is equivalent to the 2-species model, which 
is described by the reaction scheme A + B — > 2A or 
A + B — > 2B (both channels are equally probable) . The 
subsequent cyclic A^-species generalization is also called 
the A^-color cyclic voter model |Tq] . 

The rest of this paper is organized as follows. In sec- 
tion II, we discuss predator-prey dynamics in infinite di- 
mensions (mean-field) and a relationship to Nambu me- 
chanics. Section III examines interface dynamics in one 
dimension. We analyze the corresponding rate equations 
and show that spatial organization into an alternating 
mosaic of growing domains occurs for N < 5 only. While 
the qualitative predictions made in Sec. Ill are correct, 
the quantitative predictions fail. In Sec. IV, we further 
analyze the interface dynamics using primarily scaling 
arguments and numerical simulations for the most inter- 
esting cases, N = 3 and 4. We consider both sequential 
and parallel dynamics evolution rules, as the system may 
be sensitive to such rules. Section V studies the dynam- 
ics of mutations and quantities such as the fraction of 
persistent sites. In sections III-V we focus on symmet- 
ric and uncorrelated initial conditions where all the N 
species have an initial density of l/N. In section VI, we 
describe several natural generalizations of the model to 
asymmetric initial concentrations, symmetric interaction 
rules and reaction-diffusion descriptions. Section VII dis- 
cusses our results within the general framework of coars- 
ening phenomena. A Summary is presented in section 
VIII. 
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II. DYNAMICS IN INFINITE DIMENSION 

Since we are interested in the role of spatial correla- 
tions, we study primarily the extreme case of one dimen- 
sion, where spatial inhomogeneities are most pronounced. 
As a preliminary step, however, it proves useful to ex- 
amine the opposite extreme where no correlations are 
present, i.e., the cyclic Lotka- Volterra system on a com- 
plete graph. On this structure, all sites are neighbors and 
the spatial structure is irrelevant, i.e., d — oo. 

Consider the 3-species Lotka- Volterra model ([!]). The 
simplest mean-field description of this model, which is 
exact only on a complete graph, consists of the following 
rate equations for the species concentrations a(t), b(t), 
and c(t), 



a = a(b — c), 
b = b(c — a), 
c = c(a — b), 



(2) 



where the over-dot denotes the time derivative. These 
equations obey the trivial conservation law, a + b + c = 1 , 
that merely reflects particle conservation. There is an 
additional hidden integral, abc = const. The existence 
of these two integrals allows us to describe the behavior 
of the system without actually solving Eqs. @. Clearly, 
a generic solution to Eqs. (g) is periodic; this solution 
can be expressed through elliptic functions. There are 
also simpler exceptional solutions, corresponding to sit- 
uations when one of the concentrations is zero. In this 
case, the solution approaches an attracting point; e.g., if 
c = 0, one readily finds that a(t) — > 1 and bit) — > when 
t — > oo, while the approach towards these limiting values 
is exponential in time. 

Several authors notified that some Lotka- Volterra 
equations admit a Hamiltonian structure pj. However, 
even in very simple situations the Hamiltonian structure 
can be impossible, e.g., it is clearly the case when the 
number of species is odd as in the 3-species cyclic Lotka- 
Volterra system. We now describe a more intricate rela- 
tionship between the 3-species Lotka- Volterra equations 
and Nambu mechanics. Nambu mechanics |]l8| , |l9f is a 
generalized Hamiltonian mechanics where the usual pair 
of canonical variables of the Hamiltonian are replaced 
by a triplet of coordinates. Nambu formulated these dy- 
namics via the ternary operation, the Nambu bracket, re- 
placing the binary Poisson bracket. The Nambu bracket, 
{f,g,h}, satisfies the usual skew-symmetry, the Leib- 
niz rule, and the generalized Jacobi identity Let 
c\,C2,cz denote the coordinates. The canonical Nambu 
bracket is 



{/i> Si-, J3} — 



d(ci,c 2 ,c 3 ) ' 



(3) 



where the right-hand side is the Jacobian. Similarly, 
one can define the canonical V-ary Nambu bracket, 



■■■,/«) 
d(c!,...,c N ) ' 



bracket for V — 2. To define the Hamiltonian dynam- 
ics one needs a special function H; then, an arbitrary 
function f evolves according to Hamilton's equations of 
motion, / = {f,H}. To define the Nambu dynamics, 
one needs two "Hamiltonians" denoted by Hi and H 2 . 
Nambu's equations of motion read 



f = {f,H u H 2 }. 



(4) 



A few systems possessing hidden integrals of motion 
were immersed in the framework of Nambu mechanics 
p8| Surprisingly, Eqs. (Q) provide a very simple ad- 
ditional example. Indeed, writing a = c%, b = C2, c = C3, 
using the integrals of Eqs. (||) as the Hamiltonians, 
Hi = a + b + c and H 2 = abc, choosing the canonical 
Nambu bracket @), and specializing Nambu equations 
(0) to / = a, b, c, we recover Eqs. (||). 

One can investigate Lotka- Volterra models for a 
greater number of species. We checked that for N = 4, 
one can still recast the corresponding Lotka- Volterra 
equations for the concentrations Cj,l < j < 4, into the 
framework of Nambu mechanics with canonical Nambu 
bracket and the Hamiltonians Hi — ci + c 2 + C3 + C4, 
H 2 = C1C3, and H3 — c 2 c^. It would be interesting 
to establish a Nambu-type formulation (if any) for the 
cyclic V-species model with arbitrary N. On the com- 
plete graph, the corresponding rate equations are 



Cfc — c k (c k+ i — C k _i) 



(5) 



for k — 1, . . . ,N, and the addition and subtraction of 
the index are taken modulo N. Eqs. (^) appear not 
only in mathematical biology but also in plasma 
physics where they describe the fine structure of the spec- 
trum of Langmuir oscillations (22| (the system (|5|) with 
— 00 < k < 00 is also called "Langmuir lattice"). 

The possible existence of the Nambu formulation im- 
plies that Eqs. (0) should satisfy the Liouville condi- 
tion pjj, J2k dFk/dck = 0. An elementary computation 
shows that this necessary condition is indeed satisfied. 
Another necessary condition is the existence of N — 1 in- 
dependent integrals of motion. For the general V-species 
cyclic Lotka- Volterra model one easily finds two such in- 
tegrals, 



Hi = } cj , and H 2 = , 



(6) 



When the number of species is even, N = 2P, the prod- 
ucts of even and odd concentrations separately, 



p p 
H3 = Y[c 2j , H i = \\c 2o _i, 

3=1 3=1 



(7) 



which is identical to the canonical Poisson 



are also conserved so there are at least three independent 
integrals of motion: Hi, H%, and H4. Generally, [-/V/2] 
integrals of Eqs. (||) involving polynomials of Cj and lncj 
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are known (see e.g. [£4|,£5| and references therein). How- 
ever, for the iV-species cyclic Lotka-Volterra model to 
admit the Nambu formulation, N — 1 independent inte- 
grals are necessary, and it is not clear whether a Nambu 
structure exists for N > 5. 

Although we have not yet derived new results with the 
help of the Nambu structure, its richness may be useful 
in understanding the evolution of cyclic Lotka-Volterra 
systems. 



The hoping rate was taken as unity without loss of gen- 
erality. Solving this equation subject to the initial con- 
ditions 1(0) = 1/2 gives 

lit) = (9) 
w 2 + 4< w 

The system evolves into a mosaic of alternating domains 
AABBBAAABBB. The typical size of a domain grows 
linearly with time, C(t) ~ t. 



III. RATE EQUATIONS FOR ID INTERFACE 
DYNAMICS 

The infinite dimensional analysis fails to describe the 
dynamics of the actual stochastic process in low dimen- 
sions. For example, while the sum H\ = ■ Cj is con- 
served, the product H2 = ^ijCj is not a conserved quan- 
tity in ID. Furthermore, the structure of Eqs. (||) does not 
address fluctuations in the spatial distribution of inter- 
acting populations. For N < 5, we shall see that the spa- 
tial structure evolves forever, single-species domains arise 
and grow indefinitely, and the process exhibits coarsen- 
ing. In other words, equilibrium is never achieved and 
instead a networks of domains develops. The domain 
patterns are self-similar, i.e., the structure at later times 
and at earlier times differ only by a global change of scale. 
Such a behavior is a signature of dynamical scaling. 

In the following section, we study the motion of "do- 
main walls" , namely, interfaces separating domains of dif- 
ferent species. For the N species process a bond connect- 
ing two sites is an interface bond if the corresponding 
two sites are occupied by two different species. Thus, 
there are N — 1 independent types of interfaces, of which 
N — 3 are immobile and 2 are mobile. For symmetric 
initial conditions (ci(0) = 1/-/V), the different types of in- 
terface bonds are present with initial concentration equal 
to l/N (with probability 1/N a given bond does not con- 
tain an interface). Interfaces move and react according 
to TV-dependent rules defined below. For large N, most 
interfaces are immobile and the system quickly reaches a 
state where all mobile interfaces are eliminated. 



B. 3-species 

In the case N = 3, there are two types of interfaces, 
right moving (AB, BC, and CA), and left moving (BA, 
CB, and AC), denoted by R and L, respectively. Start- 
ing with a symmetric initial distribution, all right (left) 
interfaces are equivalent. When a right moving interface 
meets a left moving one, they annihilate R + L — > 0. 
When a right moving interface overtakes another right 
moving interface, they give rise to a left moving interface, 
R+R — > L, and similarly, L+L — > R. The corresponding 
rate equations are 

R = -2R 2 - 2RL + L 2 , 

L = -2L 2 -2RL + R 2 . (10) 

The interface concentration is readily found, 

m=Ht) = ^. (id 

The behavior is similar to the case N = 2, as the re- 
sulting spatial pattern form a mosaic of single-species 
domains whose average size is C ~ t. The previous anal- 
ysis implicitly assumes that interfaces hop one at a time, 
namely sequential dynamics. Alternatively, one can con- 
sider simultaneous hoping, or parallel dynamics. Here 
interfaces move ballistically, and thus, interfaces mov- 
ing with the same velocity do not interact. The reac- 
tion scheme is R + L — ► 0, and the rate equations read 
R = L = —2RL. The resulting interface concentrations 
R{t) = L(t) = l/(3+2i) differ only slightly from Eq. (0). 



A. 2-Species 

Consider the simplest case of N = 2, where there are 
two equivalent interfaces (AB and BA), denoted by /. 
An isolated interface performs a random walk, i.e., it 
hops to one of its nearest neighbors. When two inter- 
faces meet they annihilate. The corresponding reaction 
scheme is therefore / + I — > 0. Assuming that neighbor- 
ing interfaces are uncorrelated, we arrive at the binary 
reaction equation 

i = -U 2 . (8) 



C. 4-species 

In the 4-species model there appear static interfaces 
denoted by S (AC, BD, CA, and DB), in addition to 
the previously defined right moving interfaces (AB, BC, 
CD, and DA), and left moving interfaces (BA, CB, DC, 
and AD). Interfaces react upon collision according to the 
rules R + L^%, R + S^L, R + R^S,L + L^S, 
and S + L — > R, resulting in the following rate equations 

R= — 2R 2 — 2RL — RS + SL, 

L= -2L 2 -2RL - SL + RS, (12) 

S= R 2 + L 2 - RS - SL. 
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Solving these equations subject to the appropriate initial 
conditions gives 



R(t) = L(t) 



A + At' 



S(t) 



1 



y/A + H 4 + At 



(13) 



Different rules govern the decay of static and mo- 
bile interfaces, and consequently, the coarsening pro- 
cess is characterized by two intrinsic length scales. 
The typical distance between two static interfaces, i 1 / 2 , 
grows slower than the distance between two moving 
interfaces, t. A nontrivial spatial organization oc- 
curs in which large "supcrdomains" contain many do- 
mains of alternating noninteracting (AC or BD) species, 
BAAACCAAACCCAAAD. We denote the typical do- 
main size by I ~ i 1 / 2 and the typical superdomain size by 
C <~ t. The typical number of noninteracting domains in- 
side a superdomain grows as £/£ <~ t 1 / 2 . Such an organi- 
zation is a consequence of the existence of noninteracting 
species, which first occurs at N = 4. 

It is useful to consider parallel dynamics as well. 
Again, the reaction scheme is altered only in that the 
interfaces moving in the same direction do not interact. 
The reaction scheme, R + L^$, R + S^L, and 
S + L — > R, is described by the following rate equations 



R= - 2RL - RS + SL, 
L= -2RL- SL + RS, 
S = -RS-SL. 

Solving the rate equations we arrive at 

1 



(14) 



R(t) = L(t) = S(t) 



A + 2t 



(15) 



Interestingly, when N — A coarsening is sensitive to the 
details of the dynamics. Parallel dynamics is governed 
by a single length scale, in contrast with the two scales 
underlying sequential dynamics. 



D. 5-species 

In the 5-species case there are two types of sta- 
tionary interfaces, Sr (AC, BD, CE, DA, EB) and S L 
(AD, BE, CA, DB, EC), in addition to the right and 
left moving interfaces, R (AB, BC, CD, DE, EA) and L 
(BA, CB, DC, AD, AE). The reaction process is sym- 
bolized by R + L -» 0, R + S L -» L, R + S R -> S L , 
S R +L -> R, S L +L -> S R , R+R -» S R , and L+L -» S L . 
In other words, when a moving interface hits a stationary 
interface of the same kind, the outcome is a stationary 
interface of the opposite kind; when a moving interface 
hits a stationary interface of the opposite kind, a dissim- 
ilar moving interface emerges. Collisions between simi- 
lar moving interfaces produce stationary interfaces of the 
same kind, and thus, the obvious notations Sl and Sr. 
For the 5-species model with sequential dynamics, the 
rate equations read 



L = 
R = 
S L = 
Sr = 



- 2RL - LS L - LS R + RS L - 2L 2 , 

— 2RL — RSl — RSr + LSr — 2R , 
-RSl-LSl + RSr + L 2 , 
-RSr-LSr + LSl + R 2 - 



(16) 



The reaction scheme and consequently the rate equa- 
tions arc invariant under the duality transformation 
(R,S R ) < — ► (L,S L ). Particularly, for R(0) = L(0) and 
Sr(0) — Sl(0), the corresponding densities remain equal 
forever. This condition is certainly satisfied for the sym- 
metric initial conditions R(0) = L(0) = S L (0) = S R (0) = 
1/5. Therefore, R(t) = L(t) and S L (t) = S R (t), and in 
what follows we shall use the notations M(= R = L) 
for mobile interfaces and S(= Sl = Sr) for stationary 
interfaces. Thus, the four rate equations reduce to a pair 
of rate equations: 



M = -4M 2 - SM, S = M 2 - SM. 



(17) 



These equations can be linearized by introducing a modi- 
fied time variable, T(t) — J Q M(t')dt' . Using the notation 
' = d/dT, we rewrite the governing equations as 

M' = -AM -S, S' = M- S. (18) 

Solving these equations gives 

M(T) = \ (A+e-^ A + T - A_e-^ A - T ) , 
S(T) = \ (\+e-^ x - T - A_e-^ A + T ) , (19) 

with the shorthand notations, A± = (V5 ± l)/2. Here, 
the moving interfaces are depleted at = 2(lnA+)/\/5. 
The density of static interfaces approaches a finite value, 

S(oo) = i (A 2 ^ ^ - X 2 _ +VE ^ ^ 0.152477, so the average 

domain size in the frustrated state is £(oo) = l/2S'(oo) = 
3.27918. In terms of the actual time t, the density of 
moving interfaces decays exponentially, R(t) oc g-sC 00 )*. 
Contrary to the previous cases, N < 5, no coarsening 
occurs and the system quickly approaches a frozen state 
of short noninteracting same-species domains separated 
by stationary interfaces. 

A similar picture is found for parallel dynamics as well. 
Here, the reaction process is R + L — > 0, R + Sl — > L, 
R + Sr —> Sl, Sr + L —> R, Sl + L —> S R , and the rate 
equations are 



L 
R 
S L 
Sr 



2RL — LSl — LSr + RSl, 
2RL — RSl ~ RSr + LSr, 
RSl — LSl + RSr, 
RSr. — LSr + LSl- 



(20) 



The useful duality relation, (R,Sr) < — ► (L,Sl), still 
applies, so there are only two independent interface con- 
centrations, M and S, which evolve according to the fol- 
lowing rate equations 
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M = -2M 2 - MS, S = -MS. (21) 

The calculation is very similar to the sequential case, and 
we merely quote the results 

M(T)= , 5(T) = — . (22) 

The limit t — > oo corresponds to T — * — In 2. We find 
that the density of static interfaces saturates at a finite 
value, 5(oo) = 1/10, while the density of moving inter- 
faces decays exponentially in time, M(t) oc e - " 5 '- 00 ^ = 
e -t/io_ rp ne avera g e s j ze f a d oniam i n the frozen state 

is £(oo) = l/25(oo) = 5. 

To summarize, both for parallel and sequential dynam- 
ics the rate equations predict coarsening when the num- 
ber of species is sufficiently small, N < 5, and fixation for 
a large number of species, N > 5. When fixation occurs, 
each site attains a final state while for TV < 4 the state 
of any site continues to change, although the frequency 
of changes decreases with time. It is remarkable that 
the rate equation approach which neglects spatial corre- 
lations between interfaces correctly predicts the marginal 
food chain length for fixation, N c — 5, as has been proved 
rigorously for both sequential dynamics Jl6| and parallel 
dynamics |p(f . 

In the coarsening cases, N < 5, both for the two- and 
three-species case the average domain size £(t) grows lin- 
early with time independent of the dynamics. In the 4- 
species case, however, the rate equation theory predicts 
linear growth £(t) ~ t for parallel dynamics, and slower 
"diffusive" growth £(t) ~ \ft for sequential dynamics. In 
the latter case, the larger linear scale still exists and it 
characterizes the typical distance between two mobile in- 
terfaces. 



IV. COARSENING DYNAMICS IN ONE 
DIMENSION 

The above rate equation theory successfully predicts 
the fixation transition at N c = 5 in agreement with rig- 
orously known results Jll|^6| . If one assumes that the av- 
erage concentration of each species is conserved through- 
out the process, which is clearly correct at least in the 
case of equal initial concentrations, one can find a sim- 
ple argument for a lower bound on the marginal food 
chain length N c . A frozen chain consists of alternating 
domains of noninteracting species. For N = 2, 3 such a 
chain is impossible since all species interact. For N = 4, 
frozen chains are filled by either A and C species or B 
and D species thereby violating the conservation of the 
densities. [Note, however, that for N — 4 in finite sys- 
tems, density fluctuations could drive the system towards 
a final frozen configuration]. For N > 5, a frozen chain 
conserving the densities is possible and thus N c > 5. 
Given the kinetics predicted by the mean-field rate equa- 
tion approach usually proceeds with a faster rate than 



actual kinetics Iq], one can anticipate that the thresh- 
old number of different species predicted by the mean- 
field theory provides an upper bound for the actual N c , 
N c < 5. This is combined with the lower bound, N c > 5, 
to yield N c — 5. 

For N < 5, coarsening occurs and it is quite possible 
that the system develops significant spatial correlations. 
In such a case, the quantitative predictions of the rate 
equation theory are inaccurate. 

The cyclic iV-species Lotka-Volterra model is imple- 
mented in the following way. We consider a one- 
dimensional lattice of size A/" with periodic boundary con- 
ditions. Each site i of the lattice is in a given state Ni 
with Ni = A, B,C, . . . In sequential dynamics, we choose 
randomly a site and then one of its two nearest neighbors. 
If the neighbor is a predator of the chosen site, the state 
of the latter changes to the state of the predator. Oth- 
erwise, the state of the site remains the same. Time is 
incremented by 1/7V after each step. For parallel dynam- 
ics, all sites are updated simultaneously and change their 
state if one of their nearest neighbors is their predator. 
This cellular automata rule has been used in p6[ and it 
should be noted that the dynamics is fully deterministic. 
Coarsening behavior of the system depends on spatial 
fluctuations present in the initial state. For both types 
of dynamics, efficient algorithms keeping trace only of 
moving interfaces have been implemented. 

Below, we present numerical findings accompanied by 
heuristic arguments for the coarsening dynamics in one 
dimension. Again, we restrict ourselves to the symmetric 
initial concentration. In this case the average concen- 
tration of each species remains 1/N, despite the noncon- 
serving microscopic evolution rules. 



A. 2-species 

As mentioned previously, for N = 2, interfaces per- 
form a random walk and annihilate upon collision. This 
exactly soluble voter model |Q is equivalent to the one- 
dimensional Glauber-Ising model at zero tem per ature 
p7| , p8| . The interface concentration is given by ]2j| 

I(t)=e- 4t lI m+h(U)]/2, (23) 

with I n (x) the n th order modified Bessel function. In 
the limit t — > correlations are absent and the inter- 
face density lit) = 1/2 — t agrees with the short time 
limit of Eq. (||) . Asymptotically, the coarsening is much 
slower in comparison with the rate equation predictions, 
I {t) ~ (87rt) -1 / 2 . The system separates into single 
species domains as follows 

AAABBBB AAAABBAAA. (24) 

i 

The average domain size i exhibits a diffusive growth law 
lit) - t a with a =1/2. 



5 



In the 2-species case, parallel dynamics is generally 
meaningless since it quickly leads to a chain consisting 
of alternate clusters with no more than two consecutive 
similar-species sites. All sites of the chain then change 
their states at each time step. However, a parallel real- 
ization of the sequential dynamics is possible if all initial 
distances between interfaces are even integers. In this 
case, the number of interfaces gradually decreases, and a 
behavior very similar to the sequential case emerges [p9[ . 



B. 3-species 

It is convenient to consider first the simpler parallel dy- 
namics where interfaces move ballistically with velocity 
±1, and annihilate upon collisions. In this well under- 
stood ballistic annihilation process [ p6|j3C| -|3^] , a simple 
combinatorial calculation (see subsection V.C) yields the 
following interface density 



R(t) = 



1 



3 2t+l 



E 

t-1 

+E 



2t 
2i + l 



2t-2i 
t - i 



2t-2i-\ 
t-i 



(25) 



In the long time limit, the interfaces concentration de- 
cay R ~ (67rf) -1 / 2 is much slower than the t~ x decay 
suggested by the rate equation (|ll|) . The decay law gov- 
erning the interface density can be simply understood. 
Consider a finite interval of size L containing interfaces 
with initial concentration cq. The total number of inter- 
faces is N — - cqL. If the initial conditions are random, the 
difference between the number of left and right moving 
interfaces is roughly AN = \Nr — Nl\ ~ */N. At long 
times, all minority interfaces are eliminated and thus, the 
interface concentration approaches AN/L ~ (cq/L) 1 / 2 . 
By identifying the box size with the appropriate ballistic 
length L ~ Vt, the time dependent interface concentra- 
tion for an infinite system is found, R(t) ~ (co/Vt) 1 / 2 . 

The system organizes into large ballistically growing 
superdomains. Each superdomain contains interfaces 
moving in the same direction, neighboring superdomains 
contain interfaces moving in the opposite direction, etc. 
In addition to the average size of superdomains, there is 
an additional length scale in the problem corresponding 
to the distance between two adjacent similar velocity in- 
terfaces. We define these relevant length scales using the 
following illustrative configuration 



B AABBB CCCC AAABBCCC B, 



(26) 



The corresponding coarsening exponents, a and (3, are 
defined via I ~ t a and C ~ t 13 , respectively. For the 
N = 3 with parallel dynamics we thus find a = 1/2 and 
[3 = 1 . Starting from an initially homogeneous state, the 



system develops a unique spatially organized state which 
is a mosaic of mosaics. Indefinitely growing superdo- 
mains contain a growing number V = £/£ ~ t 1 / 2 of cycli- 
cally arranged domains {ABC ABC and CBACBA). 

We now turn to the complementary case of sequen- 
tial dynamics. Interfaces perform a biased random walk 
and thus, the ballistic motion is now supplemented by 
superimposed diffusion. In addition, two parallel mov- 
ing interfaces can annihilate and give birth to an oppo- 
site moving interface. It proves useful to consider the 
continuum version of the model where interfaces move 
with velocity +V and —V with equal probabilities, and 
have a diffusivity D. To establish the long-time behav- 
ior we assume that the system organizes into domains 
of right and left moving interfaces. Inside a domain, in- 
terfaces moving in the same direction can now annihi- 
late via a diffusive mechanism, unlike the parallel case. 
On slower than ballistic scales, the problem reduces to 
diffusive annihilation, X + X — > 0, where X is either 
R or L, with a density decaying as Cdiff(i) ~ (Dt)^ 1 / 2 . 
On ballistic scales the problem (almost) reduces to the 
ballistic annihilation process R + L — > 0, with the den- 
sity decay Cb a ii(£) ~ (co/Vt) 1 / 2 as described previously. 
However, to describe the complete ballistic-diffusion an- 
nihilation, one cannot use the initial concentration cq 
since it is constantly reduced by diffusive annihilation. 
Therefore, we replace the initial concentration Co with 
the time dependent concentration Cdis{t), and we find 
fS5| , [36"f c ~ (DV 2 ^)^ 1 / 4 which in particular implies 



+3/4 



(27) 



This result is quite striking since separately both an- 
nihilation processes, diffusion-controlled and ballistic- 
controlled, give the same coarsening exponent 1/2, so 
one expects that their combination does not change the 
behavior while in fact it enhances the coarsening expo- 
nent to 3/4. The resulting spatial structure is similar to 
the parallel case, Eq. (p6|). However, the smaller length 
scale is now a geometric average of a diffusive and a bal- 
listic scale as follows from Eq. ( p7| ) , while the larger scale 
remains unchanged, C ~ t. 

We have performed Monte Carlo simulations for 100 re- 
alizations on a lattice of size 10 6 , for times up to t ~ 10 6 . 
Results are shown in the Fig. 1. The interface concen- 
tration decays algebraically, R cx t~ a with an exponent 
a = 0.79. A careful analysis shows that the local slope 
a(t) = d\nc(t)/d\nt should asymptotically approach the 
value —3/4. It is possible that this finite time effect can 
be attributed to the recombination reaction (R+R — ► L), 
which is not an annihilation reaction. Nevertheless, a 
single L-interface inside an i?-domain is quickly annihi- 
lated by nearest i?-interface, and therefore recombina- 
tion is asymptotically equivalent to annihilation. We ex- 
pect that as t — > oo, the coarsening exponent is indeed 
a = 3/4. 
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FIG. 1. The concentrations of interfaces as a function of 
time for the 3-species model and a sequential dynamics in a 
log-log plot. A line of slope 3/4 is shown as a reference. The 
insert shows the local exponent a(t) as a function of 1/lni. 
A limiting value of a — + 3/4 is plausible. 

To summarize, the resulting spatial patterns in the 
N = 3 case consist of superdomains of cyclically arranged 
domains as in Eq. (||^). The larger length scale is bal- 
listic, C ~ t, while the smaller length scale is sensitive 
to the microscopic details of the dynamics: t ~ W 2 for 
parallel dynamics and I ~ i 3 / 4 for sequential dynamics. 



C. 4-species 

For the 4-species cyclic Lotka-Volterra model, numer- 
ical simulations indicate that parallel and sequential dy- 
namics are asymptotically equivalent and that the do- 
main structure is qualitatively similar to the predictions 
of the sequential rate equations, M(t) <C S(t). We use 
heuristic arguments to obtain the values of the coarsen- 
ing exponents a and j3, characterizing the density decay 
of mobile M(t) ~ t~ a , and static S(t) ~ interfaces. 

Given the equivalence of parallel and sequential dy- 
namics, we restrict ourselves to the simpler former dy- 
namics. What is the spatial structure in the long time 
limit? Since M(t) < S(t) (here and below M(t) = 
R(t) = L(t) denotes the density of moving interfaces), 
we assume an alternating spatial structure of "empty" 
regions (with no more than one moving interface) and 
"stationary" regions (with many stationary interfaces in- 
side any such region). If the interface densities obey scal- 
ing, then the size of the empty and the stationary regions 
should be comparable. The typical size of an empty or a 
stationary region is therefore of the order of M _1 . The 
typical number of stationary interfaces inside a stationary 
region is of the order of S/M. The evolution proceeds as 
follows: A moving interface hits the least stationary par- 
ticle and bounces back (since R+S — ► L and S + L — ► R). 



Then this interface hits the least stationary particle of the 
neighboring stationary region, and bounces back again. 
This "zig-zag" process continues and at some time one 
of these stationary regions "melts" , thereby giving birth 
to a larger empty region. If there is a moving particle 
inside merging empty region, the two moving particles 
quickly annihilate. If there is no such particle, the mov- 
ing particle continues to eliminate stationary interfaces. 
This process is illustrated in Fig. 2. 





FIG. 2. Space time diagram of the interface motion in the 
4-species case with parallel dynamics. 

The typical time r for a stationary region to melt is 
r = M~ x x S/M = S/M 2 . This melting time r is also 
the typical time for annihilation of a moving interface 
and thus, 



M 



M 



M 3 



(28) 



Substituting S(t) cx t~ a and M(t) cx t~ p into Eq. (§ 
we get the exponent relation 



2/3- a = 1. 



(29) 



In the next section, we introduce the mutation distri- 
bution, and find an equivalence between the fraction of 
persistent sites and the static interface density. Using 
this relation and a simple solvable example, we find the 
exponent relation 



a + (3 = l. 



(30) 



The two exponent relations therefore imply the values 
a = 1/3 and (3 = 2/3. We have simulated 100 sys- 
tems of size 10 6 up to times t ~ 10 6 . The results are 
shown in Figs. 3 and 4. We have found S(t) cx t -0 34 , 
R(t) = L(t) cx t~ - 69 for parallel dynamics, and S(t) cx 
i~ - 35 , R(t) = L(t) cx t~ - 70 for sequential dynamics. We 
conclude that the simulation results support the above 
predictions. 

As in the three-species case there are two relevant 
growing length scales. The system organizes into do- 
mains of alternating noninteracting species with a typi- 
cal size I ~ i 1 / 3 . On the other hand, active interfaces are 
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separated by a typical distance of C ~ t 2 l z , according to 
the following illustration 

c 

B 'AACCCAAA CCCC AACCAAACCC D . (31) 

t 

In any finite lattice, density fluctuations drive the sys- 
tem towards a final frozen or "poisoned" configuration, 
i.e. configuration filled by either A and C, or B and D. 
This poisoning happens when the size of superdomains 
becomes of the order of lattice size. The poisoning time 
is therefore proportional to A/" 3 / 2 for an TV-site chain. 



We stress that a 3- velocity ballistic annihilation model, 
R + L -> 0, R + S -> 0, S + L -> 0, has been recently 
investigated |3^-|3^|, and the symmetric case, R = L, has 
been solved exactly . For the special initial condition, 
i?(0) = L(Q) = 35(0)/2, a surprisingly similar behavior 
R(t) ~ S(t) <~ i~ 2 / 3 occurs. It would be interesting to 
establish a relationship between this solvable ballistic an- 
nihilation model and the interface motion in the 4-species 
process. 



D. 5-species 



2 




FIG. 3. The concentrations of the moving interfaces as a 
function of time for the 4-species cyclic Lotka-Volterra model 
with sequential dynamics (diamonds) and with parallel dy- 
namics (circles). The slope give the exponents /3 scq = 0.70 
and /3 par = 0.69. A line of slope 2/3 is shown as a reference. 



10 



o 

° ° o 
o 

° 

o "o 
°o 












°o o \ 






"" ' °%> 



FIG. 4. The concentrations of the stationary interfaces 
as a function of time for the 4-species model for sequential 
dynamics (diamonds) and parallel dynamics (circles). The 
slope give the exponents a scq = 0.35 and a par = 0.34. A line 
of slope 1/3 is shown as a reference. 



For the 5-species cyclic Lotka-Volterra model, it is well 
known that the system approaches a frozen state Jl6|,^6| . 
The kinetics approach towards saturation has not been 
established, though. 

We now present a heuristic argument for estimating the 
concentration decay of the mobile interfaces. Since the 
density of mobile interfaces rapidly decreases while the 
density of stationary interfaces remains finite we can ig- 
nore collisions between mobile interfaces. Thus we should 
estimate the survival probability of a mobile interface 
in a sea of stationary ones. There are two reactions in 
which moving interfaces survive although they change 
their type, R + Sl — * L and L + Sn.—f R. Thus, a right 
moving interface is long lived in the following environ- 
ment 

• • ■ SrSrSrSrMSlSlSlSl ■ ■ ■ (32) 

Clearly, in such configurations the zig-zag reaction pro- 
cess takes place. The moving interface travels to the right 
during a time to ~ 1/cofo, eliminates a stationary inter- 
face and travels to the left a time of order 2to, eliminates 
an interface and travels back to the right, etc. Thus, 
to eliminate A^ s interfaces, the moving interface should 
spend a time of order t ~ to X^=i * = toN s (N s + l)/2. 
Therefore, the number of stationary interfaces N s (t) 
eliminated by a moving interface scales with time as 
N s (t) ~ V co vat. Configurations of the type ( |32| ) are 
encountered with probability cx e~ Ns with N s the config- 
uration length, and thus, the density of moving interfaces 
exhibits a stretched exponential decay, 

M ( ( ) (xe -co. s t.x v W i (33) 

The stretched exponential behavior ( |33| ) is expected to 
appear for arbitrary N > 5. When the number of inter- 
faces exceeds the threshold value, N > 5, stationary in- 
terfaces of "intermediate" types arise, the crossover from 
initial exponential behavior to the asymptotic stretched 
exponential behavior is shifted to larger times and there- 
fore harder to observe numerically. For the threshold 
number of species, however, we have found a convincing 
agreement between the theoretical prediction of Eq. ( |3^ ) 
and numerical results (see Fig. 5). Finally we note that 
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the actual kinetics ( |33| ) is slower than the mean-field 
counterpart, Mmft(^) oc e - *, due to spatial correlations. 
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FIG. 5. The concentration of moving interfaces as a func- 
tion of the square root of time in linear-log plot for the 
5-species cyclic Lotka-Volterra model with sequential dynam- 



V. DYNAMICS OF MUTATIONS 

Consider a lattice site occupied by some species, say 
A. What is the probability that this site has been occu- 
pied by the same species during the time interval (0,i)? 
Otherwise, what is the fraction of A sites which never 
"mutated"? We denote the fraction of "persistent" A 
species by Ao(t); Bo(t), Co(t), etc. are defined analo- 
gously. We can further generalize these probabilities to 
define, e.g., A n (t), the fraction of sites that have under- 
gone exactly n mutations during the time interval (0,i). 
We start by analyzing these quantities on the mean-field 
level and then describe exact, scaling, and numerical re- 
sults in one dimension. 



A. Mean-Field Theory 

Let us investigate the 3-species cyclic Lotka-Volterra 
on the complete graph; the generalization to the N- 
species case is straightforward. The rate equations de- 
scribing the evolution of the mutation distribution A n (t) 
are 



A 3n = aA 3 „_i - cA 3n , 

A^n+l — cA^n — bA^ n+ \, 
A?, n +2 = bA^n+i — aAy, n+ 2 



(34) 



with A-\(t) = 0. Analogous equations can be writ- 
ten for B„(t) and C n (t) by cyclic permutations. These 



rate equations form an infinite set of linear equations 
with a(t),b(t), and c(t) as (time-dependent) coefficients. 
Therefore, the general case is hardly tractable analyti- 
cally since the coefficients, i.e. solutions of Eqs. (0), are 
elliptic functions. We therefore restrict our attention to 
the symmetric case a = b = c = 1/3, and examine P n (t), 
the total fraction of sites mutated exactly n times. The 
quantity P n (t) evolves according to 



Pti — Pn- 



Pn, 



(35) 



with P_i = to ensure Pq = —Po- In Eq. ( |3q ) we ab- 
sorbed the concentration factor 1/3 into the time-scale 
for convenience. Solving ( |35| ) subject to the initial con- 
dition P n (0) = 5 n o, one finds a Poissonian mutation dis- 
tribution 



Pn{t) 



t" 



(36) 



This mutation distribution is identical to the one found 
for the voter model |ll| on the mean-field level. 

The distribution is peaked around the average (n) = t, 
and the width of the distribution, a, is given by a 2 = 
(n 2 ) — (n) = t. In the limits, t — > oo, n — > oo, and 
(n — t)/y/i finite, P n (t) approaches a scaling form 



Pn(t) 



1 



(37) 



where the scaling distribution function $oo(z) is Gaus- 
sian, $00(2;) = (27r) -1 / 2 exp(— z 2 /2), and the index 00 
indicates that the solution on the complete graph cor- 
respond to the infinite-dimensional limit. We also note 
that the fraction of persistent sites decreases exponen- 
tially, P (t) = e"*. 

Let the initial state of a site be A, without loss of gen- 
erality, then the probability that the state is A at time 
t is given by Ro(t) = ^2 n P3n(t). In general, three such 
autocorrelation functions 



(38) 



n=0 



correspond to the three possible outcomes at time t, A 
if jfe = 0, C if k = 1, and B if k = 2. The quantity 
i?o(t) is evaluated from equation ([36]) using the identity 
e* + e ct + e cH = 3 J2„ t 3n /(3n)\, with C = e 27Ti / 3 . Gener- 
ally, we find that 



Rkit) = - 



l + 2e- 3t / 2 cos(^+^ 



(39) 



for k — 0, 1, 2. The structure of the autocorrelation func- 
tions is rather simple - an exponential approach to the 
equilibrium value Rk (00) = 1/3 is accompanied by oscil- 
lations. The three autocorrelation functions differ only 
by a constant phase shift. One can verify that expo- 
nential decay occurs for arbitrary N, and that temporal 
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modulations occur when N > 2. Numerical simulations 
for the case N — 3 indeed show temporal oscillations in 
Ro(t). However, algebraic rather than exponential decay 
is found for the magnitude. 



B. Scaling behavior 

Mutation dynamics and coarsening dynamics are 
closely related [^3|. For example, the rate of mutation 
is given by the density of moving interfaces. Using simi- 
lar scaling arguments, we study asymptotic properties of 
the mutation distribution in the one-dimensional case. 

The mutation distribution satisfies the normalization 
condition, ^ n P n = 1. Let the average number of mu- 
tations be (n) = Y] n nP n . Every motion of an interface 
contributes to an increase in the number of mutations in 
one site, and thus the mutation rate equals the density of 
moving interfaces, d(n(t))/dt — M(t). In the coarsening 
case, TV < 5, we found that the moving interface density 
decays algebraically, M[t) ~ t~^. Therefore, the average 
number of mutations grows algebraically, (n(t)) ~ f ', 
with v = 1 — fi. For N — 2 and 3, the density of moving 
interfaces decays inversely proportional to the average 
domain size, M ~ I , since stationary interfaces are ab- 
sent when N < 3; therefore, fj, = a. For N = 4, however, 
the density of moving interfaces is inversely proportional 
to the average size of superdomains, M ~ C , implying 
(1 = 0. 

In the case of N = 2, it has been shown that the muta- 
tion distribution obeys scaling . We assume that this 
behavior generally holds when the system coarsens, 



Pnit) 



1 



■ $ 



n(t)> \{n{t)) 



(40) 



The behavior of the scaling function, in the 

limit of small and large arguments z reflects the frac- 
tion of persistent and rapidly mutating sites, respectively. 
Typically, the fraction of persistent sites decays alge- 
braically in time, Po(t) ~ t~ e , with 6 the persistence 
exponent. This exponent has been studied recently in 
several contexts such as kinetic spin systems with con- 
servative and non-conservative dynamics and diffusion- 
reaction systems ]37|-^|. In the N — 2 case, the scaling 
function was found to be algebraic, <&(z) ~ z 1 , in the 
limit z — ► 0. Assuming this algebraic behavior for N = 3 
and 4 as well implies the exponent relation 



61 = 1/(7 + 1). 



(41) 



The large z limit describes ultra-active sites. A con- 
venient way to estimate the fraction of such sites is to 
consider sites which make of the order of one mutations 
per unit time. At time t, the fraction of these rapidly mu- 
tating sites is exponentially suppressed, Pt(t) oc exp(— t). 
It is therefore natural to assume the exponential form 
$(z) ~ cxp(— z s ) for the tail of the scaling distribution, 



thereby implying an additional exponent relation /i<5 = 1 . 
To summarize, the scaling function underlying the mu- 
tation distribution has the following limiting behaviors 



f z~< z<l; 
1 exp(— const, x z s ) z 1. 



(42) 



In section IV, we obtained the exponent /i, charac- 
terizing the decay of moving interfaces. The mutation 
exponent and the tail exponent are readily found us- 
ing the respective exponent relations, v = 1 — fi and 
5 = l//i. To find the persistence exponent, we note 
the equivalence between the fraction of persistent sites 
and the fraction of unvisited sites in the interface pic- 
ture mif]. For N = 2, the value 6 = 3/8 has been 
established analytically J39|. For N = 3, different be- 
haviors were found for parallel and sequential dynam- 
ics, and therefore it is necessary to distinguish between 
the two cases. As mentioned above, the parallel case 
reduces to a two- velocity ballistic annihilation process. 
The probability that a bond has remained uncrossed from 



f 



-1/2 



SCO 



the left by right-moving interfaces is S+(t) 
Eq. (|25|). Analogously, S-(t) ~ i^ 1 / 2 , and consequently 
P (t) = S-(t)S+(t) ~ t- 1 or 6 = 1 follows |§. In the 
sequential case, we have not been able to determine the 
persistence exponent analytically, and a preliminary nu- 
merical simulation suggests that 9 = 1 as in the parallel 
case. 

For N = 4, the number of unvisited sites is equiva- 
lent asymptotically to the survival probability of a static 
interface, Po(t) ~ S(t), and using the definitions of sec- 
tion IV, we find 6 = a, i.e. 8 = 1/3. We now present 
a heuristic argument supporting the exponent relation 
(pp|). Substituting the previously established exponent 
relations v = 1 — fi, fj, = /3, and 9 = a in Eq. (J4l|) yields 



a=(l-/3)( 7 +l). 



(43) 



We now argue that 7 = and thus Eq. (|43| ) reduces 
to a + P = 1, i.e. to Eq. (|30|). We first recall that inter- 
faces in the 4-species case react according to R + S — * L, 
L + S — * R, and R + L — > 0. In the long-time limit, the 
zig-zag reactions R + S — » L and L + S — ► R dominate 
over the annihilation reaction R + L — > 0. We, therefore, 
consider a simpler solvable case where a single mobile 
interface is placed in a regular sea of static interfaces. 
Then this interface moves one site to the right, two to 
the left, three to the right etc. In a time interval (0,i), 
this interface eliminates N s ~ i 1 / 2 static interfaces. The 
origin is visited iV s times, site 1 is visited N s — 1, site 
— 1 is visited N s — 2, etc. This implies that the mutation 
distribution is P n (t) = (n) _1 $(n/(n)), with (n) ~ N s 
and &(z) = 1 for z < 1 and = for z > 1. Hence, 
ignoring the annihilation reaction leads to 7 = 0. This 
approximation is inappropriate for predicting the tail of 
$(z) which is sensitive to annihilation of the moving in- 
terfaces. However, in the small z limit the annihilation 
process should be negligible, and thus 7 = 0. 
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Monte-Carlo simulations confirm the anticipated scal- 
ing behavior of Eq. (|40|). In Fig. 6, the scaled muta- 
tion distribution function (n)P n (t) is plotted versus the 
scaled mutation number n/{n), for a representative case 
N = 4 at different times t = 10 3 ,10 4 ,10 5 . It is seen 
that the plots are time independent. Furthermore, the 
scaling function approaches a finite nonzero value in the 
limit of small z = n/{n), in agreement with the scaling 
predictions, 7 = 0. 



0.8 




n/<n> 

FIG. 6. The mutation distribution scaling function for a 
representative case of N = 4 with sequential dynamics. Sim- 
ulations of 10 realizations of a system of size 10 6 for time 
t — 10 3 (circles), t — 10 4 (squares) and t — 10 s (diamonds). 

In summary, coarsening dynamics can be character- 
ized by a set of exponents a, P,^,8,v,6. Table 1 gives 
the values of these exponents which are believed to be 
exact, although for some of the exponents only numeri- 
cal evidence exists so far. 



N 


a 


P 


V 


5 


6 


7 


2 


1/2 




1/2 


2 


3/8 


-1/4 


3 (parallel) 


1/2 


1 


1/2 


2 


1 


1 


3 (sequential) 


3/4 


1 


1/4 


4 


1 


1/3 


4 


1/3 


2/3 


1/3 


3 


1/3 





4 (symmetric) 


3/8 


1/2 


1/2 


2 


3/8 


-1/4 



Table 1 Coarsening and mutation exponents in ID. 



C. An exactly solvable case 

The 3-species Lotka-Volterra model with parallel dy- 
namics is equivalent to the exactly solvable two-velocity 
ballistic annihilation |3(|. We exploit this equivalence 
to compute analytically the mutation distribution. A 
species in a given site mutates each time it is crossed by 
an interface. As the fraction of persistent sites is equiv- 
alent to the fraction of uncrossed bonds, the fraction of 



sites visited n times equals the fraction of bonds crossed 
exactly n times by the interfaces. In the symmetric case, 
the initial concentration of moving interfaces of velocity 
+1 or —1 is 1/3 (interfaces are initially absent with prob- 
ability 1/3). Interfaces move ballistically and the system 
is deterministic, i.e., any late configuration is a one-to- 
one function of the initial configuration. It is also nat- 
ural to consider integer times t. The distribution P n {t) 
for a given site is completely determined by the initial 
distribution of the interfaces on the t bonds to the left 
of this site and on the t bonds to the right of this site 
since, further interfaces can not reach the site in a time 
t. This 2t initial bonds can be mapped onto a random 
walk with uncorrelated steps of length ±1 or zero since 
interfaces are initially uncorrelated. We set Sq = and 
define Si recursively via Si = Sj-i+M,', i — 1, . . . , t where 
Vi = ±1 is the velocity of the i th interface to the right of 
the considered site and Vi = if the interface is absent. 
Similarly, S-i — S-u—i) — i = 1, ...,t. Thus, one 
has two random walks starting from the origin, (i, Si) 
and (— i, S-i), i = 0, . . . , £, with i being a time-like vari- 
able and Si the displacement. The crucial point is that 
the number of interfaces crossing the target site at the 
origin during the time interval (0, t) is given by the abso- 
lute value of the minimum of the combined random walk 
(i,Si),i = -t,...,t (see Fig. 7). 




FIG. 7. Mapping of the initial distribution of the species 
on a random walk. The absolute minimum reaches by the 
random walk is equal to the number of mutations undergoes 
by the center site. 

Indeed, the minimum attained by the random walker 
on the left (right) gives the excess of interfaces coming 
from the left (right) not destroyed by other left (right) 
interfaces that would cross the considered site. Thus, 
Pnit) equals the probability that the minimum of two 
independent i-steps random walks starting at So = is 
-n, 

n 

P n (t) = 2Q n {t) Qk(t) - Q n (t) 2 . (44) 

fc=0 
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The probability Q n (t) that a i-steps random walk start 
ing at the origin has a minimum at — n, is given by 44 



Qn(t) = Q n (t) + Q n+1 (t), 



with 



Quit) 



t—n 



i=0 



' t+n — i \ | / t — n — i N 
v 2 / ■ \ 2 / 



(45) 



(46) 



The trinomial coefficient in the above sum is set to zero if 
(t — n — i)/2 is not an integer. The sum in the right-hand 
side of Eq. ( f44j ) gives the probability that the other walker 
has its minimum at — fc, with k < n, the factor 2 reflects 
the fact that there are two random walkers. We subtract 
the last quantity Q n {t) 2 which has been counted twice in 
the summation. In particular, we have Po(t) = Qo(t) 2 , 
in agreement with the argument of the previous section. 
Also, one can show that the density of moving interfaces 
can be expressed via Qo as R(t) = Q n (2t)/3, leading to 
Eq. (§§. 

To determine the asymptotic behavior of P n (t) we first 
compute Q n (t) . Making use of the Gaussian approxima- 
tion for the trinomial coefficients we find 



(47) 




Combining (47) 



and (44) yields 



P„(t)~ A /-Erf 



-3n 2 / 4 * 



(48) 



with Erf(z) = J Q Z due " 2 . The existence of an exact 
solution is very useful for testing the validity of the scal- 
ing assumption. Indeed, Eq. ( f48| ) agrees with the general 
scaling form of Eq. (|40|), and the corresponding scaling 
function is 



'Erf(z), 



(49) 



with the scaling variable z — n /y/4t/3. The limit ing be- 
havior of this scaling function agrees with the predictions 
of Eq. (fH) as well, 



$(z) ~ { 



z < 1; 
z > 1. 



(50) 



The corresponding values of the scaling exponents v — 
/.t = 1/2, 5 = 2, = 1, and 7=1, are in agreement with 
Table 1. 



VI. EXTENSIONS 

The cyclic lattice Lotka-Volterra model can be gener- 
alized in a number of directions. A natural generalization 
is to higher dimensions. Two-dimensional case seems to 



be especially interesting from the point of view of math- 
ematical biology. In the exactly solvable N = 2 case 
(the voter model), coarsening occurs for d < 2 p2[ , for 
the marginal dimension d = 2, the density of interfacial 
bonds decays logarithmically, c{t) ~ 1/ lni |ll|], while for 
d > 2, no coarsening occurs and the system reaches a 
reactive steady state. In two dimensions, our numerical 
simulations indicate that there is no coarsening, i.e. the 
density of reacting interfaces saturates at a finite value. 
For sufficiently large number of species the fixation is ex- 
pected but we could not determine the threshold value, at 
least up to N = 10 we have seen no evidence for fixation. 

Below, we mention few other possible generalizations 
and outline some of their attendant consequences. 



A. Asymmetric Initial Distribution 



We consider uncorrelated initial conditions with un- 
equal species densities. Even in the 2-species situation, 
the behavior is surprisingly non-trivial. In particular, 
the densities of both species remain constant; the per- 
sistence exponent 9a decreases from 1 to as the ini- 
tial concentration ao increases from to 1 [|13||39 , with 
9 a = 9b — 3/8 for equal initial concentrations p9|. 

Turn now to the 3-species case and consider first par- 
allel dynamics. In general, the densities of right and left 
moving interfaces are equal as well. However, the initial 
interface distribution is correlated in the general asym- 
metric case and therefore the equivalence to ballistic an- 
nihilation is less useful. We find numerically that the 
interface density exhibits the same decay as in the sym- 
metric case, c(t) ~ t -1 / 2 . To illustrate this property let 
us consider the following example where the initial densi- 
ties are ao = 1 — 2e and &o = c o = e , with e — > 0. Initially, 
the A-species dominates over the two minority species. 
While isolated £?'s are immediately eaten by the neigh- 
boring A's, C-species domain arise and soon the C's dom- 
inate the system. However, the ultimate fate of the sys- 
tem is determined by pairs of nearest neighbors which are 
dissimilar minorities, i.e. BC and CB. Initially, these 
interfaces are present with probability e 2 ; clearly, they 
are long-lived right and left moving interfaces. These in- 
terfaces are uncorrelated and thus their density decays as 
t~ x l 2 . We also performed numerical simulations for the 
3-species cyclic Lotka-Volterra model with sequential dy- 
namics, and the interface decay, c(t) ~ £ -3 / 4 , was found 
similar to the symmetric case. The interface concentra- 
tion does not provide a complete picture of the spatial 
distribution. The main difference with the 3-opinions 
voter model is that the species densities are not con- 
served, and they exhibit a more interesting behavior (see 
Fig. 8). It is possible that the limit where one species 
initially occupies a vanishingly small volume fraction is 
tractable analytically, similar to recent studies Jl^,Q of 
Glauber and Kawasaki dynamics. 
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predators and they overtake the £>'s. The A's are eaten 
first but once the C's dominate the D's, A's have less 
and less predators. The concentration of D-species and 
the density of the moving interfaces decay exponentially 
and, therefore, the system quickly reaches a frozen state 
where af 10Z0 n = 1/4 and Cf roZ cn = 3/4 (see Fig. 9). These 
constants can be simply understood. Consider only the 
initial distribution of C and D. Regions between a pair 
successive C (such regions are present with probability 
1/4) will be filled by A's. Regions between a pair of D 
as well as regions between a C and a D (present initially 
with probability 3/4) will become C domains. 
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FIG. 8. The species densities as a function of hit for the 
3-species cyclic model with unequal initial densities (ao = 0.9 
and bo — Co — 0.05). 
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FIG. 9. The species densities as a function of time 
for the 4-species cyclic model with unequal initial densities 
(a = 0.997 and b = c = do = 0.001). 

Consider now the 4-species model. Numerically, we ob- 
served a rich variety of different kinetic behaviors. Rather 
than giving a complete description, we restrict ourselves 
to a few remarks based on simulation results and heuristic 
arguments. First, the species densities are not conserved 
globally, in contrast with the symmetric initial conditions 
or the ordinary 4-opinions voter model. Furthermore, if 
the initial densities are different, the system can fixate 
and thus reach a state such as AAACCCACCA where 
the evolution is frozen. In order to illustrate the rich 
behavior of this system we consider the following initial 
conditions ao = 1 — 3e and bo = cq = do = e with e — > 0. 
Eaten by the dominant A's and with almost no preys, B's 
quickly disappear from the system. The D's are grow- 
ing because they have much food and almost no preda- 
tors. After a while, the C's also have some food and no 



B. Symmetric Rule 

Let us now consider the iV-species Lotka-Volterra 
model with a symmetric eating rule, namely we assume 
that the i th species can eat species i — 1 mod N as well 
as i + 1 mod N. 

For N = 3, all different species can eat each other 
without any restriction. This model is thus equivalent to 
the 3-opinions voter model (also called the stepping stone 
model) . In one dimension, the concentration of interfaces 
is known to decay as £ -1 / 2 , see e.g. Q. 

For N = 4, the situation is more interesting since e.g. 
A can eat both B and D but cannot eat C. Thus this 
model is different from the 4-opinions voter model or the 
4-species cyclic Lotka-Volterra model. There are moving 
interfaces M between species A and B, B and C, C and 
D, and D and A, and stationary interfaces S between 
species A and C and species B and D. Each moving 
interface is performing a random walk. When a moving 
interface meets a stationary one, the latter is eliminated, 
M + S — » M; if two moving interfaces meet, they either 
produce a stationary interface M + M — > S or annihi- 
late M + M — ► according to the state of the underlying 
species. On the mean-field level, this process is described 
by the rate equations 



M = -4M 2 , S = M - SM. 



(51) 



Eqs. fl5l|), supplemented by the initial conditions M(0) 
1/2 and S(0) = 1/4, are solved to yield 



M(t) = 



At' 



12 (l + 2t)V4 3 + 6i ' 



(52) 



implying the existence of two scales, I ~ t 1 / 4 and C ~ t. 

Fortunately, an exact analysis of the 4-species Lotka- 
Volterra model with the symmetric eating rule is possible. 
Moving interfaces do not feel the stationary ones and they 
are undergoing diffusive annihilation. As a result, their 
concentration decays according to M(t) ~ i -1 / 2 . Follow- 
ing the discussion in the previous section, the fraction 
of stationary interfaces surviving from the beginning is 
proportional asymptotically to the fraction of sites which 
have not been visited by mobile interfaces up to time t, 
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S(t) ~ P (t) ~ t~ 3 / 8 §|. We should also take into 
account creation of stationary interfaces by the annihi- 
lation of moving interfaces. This process produces new 
stationary interfaces with rate of the order —dM/dt so 
the density of stationary interfaces satisfies the rate equa- 
tion 



dS 
dt 



dPo 
dt 



dM 
~dT' 



(53) 



Combining Eq. fl53j) with P (t) ~ t~ 3 ' 8 and M(t) ~ 
i -1 / 2 , we find that interfaces which survive from the be- 
ginning provide the dominant contribution while those 
created in the process M + M —> S contribute only to a 
correction of the order i -1 / 8 



5(i)-^ 3/8 h+o(t- 1/s ) 



(54) 



Thus a two-scale structure of the type ( 31 ) emerges with 
typical length , t ~ i 3 / 8 and C ~ t 1 ! 2 . The exponents 
for the 4-species Lotka-Volterra with symmetric rules are 
summarized in Table 1. These asymptotic results agree 
only qualitatively with the rate equations predictions. 
Simulation results are in an excellent agreement with 
these predictions, M(t) ~ t~ - 50 and S(t) ~ t~ - 35 (see 
Fig. 10). Refined analysis which makes use of the ex- 
pected correction of the order 0{t~ x / 8 ) enables a better 
estimate for the decay of stationary interfaces, namely 
S(t)~t~ 037 . 
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FIG. 10. The concentrations of stationary (diamonds) and 
moving (circles) interfaces as a function of time for the 
4-species model with a symmetric sequential dynamics. Lines 
of slope 1/2 and 3/8 are shown as references. The insert 
shows t 3 / a S(t) as a function of t -1 / 8 where a straight line is 
expected. 

The TV = 5 case with symmetric eating rules can be 
easily analyzed on the level of the rate equations. We 
omit the details as the analysis is similar to the one pre- 
sented in subsection III.D for the cyclic model. The con- 
clusion is similar as well, namely the system approaches 



a frozen state consisting of noninteracting domains. Ar- 
guing as in the cyclic case we conclude that the threshold 
number of species predicted by the mean-field rate equa- 
tion approach is exact, N c = 5, in agreement with our 
numerical simulations. 

We also found that the rate equation approach does 
not provide a correct description of the decay of the mo- 
bile interfaces: Mmft(£) oc e _t , while in the actual pro- 
cess M(t) oc e~* with n close to 1/4. An upper bound, 
7i < 1/3, can be established by comparing to the trapping 
process, M + T — ► T ]45| , where the survival probability 
a(t) for a particle diffusing in a sea of immobile traps, 
a{t) oc exp(— t 1 / 3 ) p5| , provides a lower bound for our 
original problem, Mjt) > u{t). 



C. Diffusion- Reaction Description: Cyclic Models 

So far we studied population dynamics occurring on a 
lattice. Although similar descriptions has been used in 
several other studies J7|-pd]|, the diffusion- reaction equa- 
tion approach is more popular [|l]j3[Q]. It is therefore 
useful to establish a relationship between the two ap- 
proaches. 

To this end, consider a 3-species system with particles 
moving diffusively and evolving according to the reac- 
tion scheme (|l|), supplemented by reproduction and self- 
regulation. On the level of a diffusion-reaction approach, 
this process is described by the following partial differ- 
ential equations 



dx 2 
d 2 b 



+ a(l — a) + ka(b — c), 



da 
di 



(55) 



dc 
dt 



8 2 c 
dx 2 



+ c(l - c) + kc(a - b). 



In these equations, a = a(x, i), b — b(x, t), and c = c(x, t) 
denote the corresponding densities at point x on the 
line; a(l — a) is the Lotka term describing reproduction 
and self-regulation; the diffusion constant and the growth 
rates of each species are set equal to unity, and the con- 
stant k measures the strength of the competition between 
species. 

For noninteracting species, k = 0, and Eqs. (55) decou- 



ple to the well-known single-species Fisher-Kolmogorov 
equations This equation has two stationary solu- 

tions, a — and a = 1; the former is unstable while the 
latter is stable so any initial distribution approaches to- 
ward it. Starting from an initial density close to stable 
equilibrium for x < and to unstable equilibrium for 
x > 0, a wave profile is formed and moves into the unsta- 
ble region j|,[3]4(|. The width of the front is finite as a 
result of the competition between diffusion which widens 
the front and nonlinearity which sharpens the front. 

Consider now the case of interacting species, k > 0. 
The initial dynamics is outside the scope of a theoretical 
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treatment and should be investigated e.g. numerically 
solely on the basis of Eqs. (fsB]). However, as the coars- 
ening proceeds, single-species domains form. Inside say 
an A-domain, the density of A species is almost at stable 
equilibrium, a(x,t) = 1, while the densities of B and C 
species are negligible. In the boundary layer between say 
an A and a B domains, the density of C species is neg- 
ligible. Domain sizes grow while the width of boundary 
layer remains finite. Therefore, in the long time limit one 
can treat boundary layers as (sharp) interfaces which are 
expected to move into "unstable" domain. 

To determine the velocity v of the interface and the 
density profiles we employ a well-known procedure |^],fl6| ■ 
Consider an interface between say an A domain to the 
left and a B domain to the right. We look for a wave-like 
solution, 

a(x,t) = a(£), 6(x,t) = &(£), £ = x-vt. (56) 

Substituting (p6|) into Eqs. ( [sB]) we arrive at a pair of or- 
dinary differential equations for the density profiles a(£) 
and &(£). These equations cannot be solved analytically. 
To determine the interface velocity, let us consider the 
densities far from the interface (£ = 0), say for £ 3> 1. 
In this region a(£) <C 1 and &(£) = 1. Therefore the 
equation for a(£) simplifies to 



+ va' + (1 + k)a = 0, 



(57) 



-AS 



into Eq. ( |57[ ) we get A — vX + 



In principle any velocity v > v m i n , with 



where a' — da/d£, etc. By inserting an exponential so- 
lution, a(£) 
(1 + k) = 0. 

Vraxn = 2\/l + k, is possible. This resembles the situation 
with the Fisher-Kolmogorov equation . According to 
the "pattern selection principle" the minimum ve- 

locity is in fact realized for most initial conditions. The 
pattern selection principle is a theorem for the Fisher- 
Kolmogorov equation (where the precise description of 
necessary initial conditions is known) [fl6f while for many 
other reaction-diffusion equations the pattern selection 
principle has been verified numerically 

Thus, for 3-species cyclic Lotka-Volterra model in ID 
we have established an asymptotic equivalence between 
the diffusion-reaction approach and the lattice one with 
the parallel dynamics. Given that the density of inter- 
faces decays as t -1 / 2 , one can anticipate the same behav- 
ior for the diffusion-reaction model. This result may be 
difficult to observe directly from numerical integration of 
the nonlinear partial differential equations (|55|), and es- 
tablishing the complete relationship between lattice and 
diffusion-reaction approaches remains a challenge for the 
future. 



D. Diffusion- Reaction Description: Symmetric 
Models 



Consider now the 3-species symmetric Lotka-Volterra 
on the level of the diffusion-reaction description. Rate 



equations of the type Eqs. ( Ba ) are useless in this case 
since they do not contain terms describing interactions 
among species. Nevertheless it proves useful to consider 
a similar symmetric system where interacting species mu- 
tually annihilate upon collision. The governing equations 
read 



d 2 a 



da 

db _ d 2 b 
di'dx 2 
Oc d^c 

- = — + c(l-c)-kc(a + b). 



— = q^2 + «(! - a ) - ka ( b + c )' 
+ 6(1 - b) - kb(c + a), 



(58) 



We again restrict ourselves to the late stages where a 
well-developed domain structure has been already formed 
[^,^8|. To simplify the analysis further we assume that 
the competition is strong, k — > oo, so neighboring do- 
mains act as absorbing boundaries. We employ a qua- 
sistatic approximation, i.e. we neglect time derivatives 
and perform a stationary analysis in a domain of fixed 
size, and then make use of those results to determine the 
(slow) motion of the interfaces. Inside say an A domain 
the density a(x) satisfies 



8 2 a 
dx 2 



+ a(l - a) =0. 



(59) 



One should solve Eq. (59) on the interval (0, L) subject 
to the boundary conditions a(0) = a(L) = 0. The size 
of the domain, L, is assumed to be large compare to the 
width of the interface, i.e., L 1. In this limit, the flux 
of A species through the interface is equal to p8[ 



F(L) 



1 



v/3 



const, x e 



(60) 



Clearly, if we have neighboring Li-domain and L2- 
domain, then the smallest of the two domains shrinks 
while the largest grows, and the interface moves with ve- 
locity F{Li) - F(L 2 ) oc e~ L2 - e~ Ll . Thus the average 
size grows according to 



-(L)cxexp(-(L)) 



(61) 



which is solved to yield (L) ~ \nt. We see that coarsen- 
ing still takes place, but it is logarithmically slow. 

Moreover, the determination of the complete domain 
size distribution can be readily performed, at least nu- 
merically. Clearly, in the late stage all sizes are large, 
L 1. Thus, only the smallest domain shrinks and 
the two neighboring domains grow while other domains 
hardly move at all. This provides an extremal algo- 
rithm: (i) The smallest domain L m i n is identified; (ii) If 
the nearest domains, L\ and L2, contain similar species, 
both interfaces are removed and a domain of length 
Li + -Lmin + L2 is formed; (iii) If the nearest domains 
contain dissimilar species, the two interfaces merge and 



15 



form a new interface at the midpoint, and thus domains 
of size Li + L m i n /2 and L2 + L m - m /2 are formed. This pro- 
cess is identical to the 3-state Potts model with extremal 
dynamics EJ . Similar one-dimensional models with ex- 
tremal dynamics have been investigated in a number of 
recent studies d^-HJ. 

Finally, we briefly discuss the symmetric rule model, 
with N = 2 and 4. In the two-species case |^,|8|, a mo- 
saic of alternating A and B domains is formed. The late 
dynamics is again extremal: the smallest domain merges 
with the two neighboring domains. The domain-size dis- 
tribution function p9| and several correlation functions 
have been computed analytically. In the 4-species 
model, there are two independent kinetic processes: Slow 
coarsening which takes place when domains of interacting 
species are neighbors, and fast mixing of neighboring do- 
mains of non- interacting species. Clearly, almost all do- 
mains will be soon populated by pairs of non-interacting 
species, so a mosaic of alternating AC and BD domains 
is formed, and the resulting dynamics of the symmetric 4- 
species model is identical to the dynamics of the 2-species 
model. 

Thus, in the symmetric case the reaction-diffusion ap- 
proach provides very different results compare to the lat- 
tice process. It was also seen that mapping onto a system 
of intervals with extremal dynamics does provide an ef- 
fective way to analyze the long time behavior. 



VII. DISCUSSION 

In this paper, we investigated one-dimensional Lotka- 
Volterra systems and found that they coarsen when the 
number of species is sufficiently small, N < 4. Typically, 
coarsening systems exhibit dynamical scaling with a sin- 
gle scale J53|. When scaling holds, analysis of the sys- 
tem is greatly simplified, e.g., the single scale grows as a 
power law, £(t) ~ t a , with the exponent a independent of 
many details of the dynamics, usually even independent 
of the spatial dimension |33|. In contrast, for the Lotka- 
Volterra models we found that the coarsening depends 
on details of the dynamics. There are two characteris- 
tic length scales: the average length of the single-species 
domains, t(t) ~ t a , and the average length of superdo- 
mains, C(t) ~ t" . Precise definition of superdomains de- 
pends on the number of species N: For N = 3 interfaces 
between neighboring domains move ballistically and su- 
perdomains are formed by strings of interfaces moving in 
the same direction; for ./V = 4, neighboring domains are 
typically noninteracting, and superdomains are separated 
by active interfaces. 

Dimensional analysis provides additional insight into 
the existence of more than one scale. Consider for sim- 
plicity parallel dynamics, where the relevant parameters 
are the initial interface concentration cq, the interface 
velocity v, and time t. There are only two independent 
length scales, Cq 1 and vt, and using dimensional analysis 



one expects 

£(t) = vtip(covt), C(t) = vt^{c Q vt). (62) 

If simple scaling holds, the length Cq 1 set by initial condi- 
tions should be irrelevant asymptotically. Thus, the scal- 
ing functions ^>{z) and ^(z) should approach constant 
values as z = covt — > 00 implying £(t) ~ C(t) ~ vt. In 
contrast, for the 3-species Lotka-Volterra model we found 
tp(z) ~ z -1 / 2 when z — ► 00. For the 4-species Lotka- 
Volterra model both scaling functions exhibit asymptotic 
behavior different from the naive scaling predictions, 
ip(z) ~ z~ 2 / 3 and &(z) ~ z -1 / 3 . For the Lotka-Volterra 
model with symmetric eating rule interfaces diffuse and 
thus the relevant length scales are Cq 1 and \J~Dt. Here, 
l{t) = s/Diip{clDt) and £(t) = \fDt^{clDt). When 
N = 4, the two scale structure implies ip(z) ~ z -1 / 8 as 
z — > 00 . 

Thus simple dynamical scaling is violated for the one- 
dimensional Lotka-Volterra models. Violations of scaling 
have been reported in a few recent studies of coarsen- 
ing in one- and two-dimensional systems p4]-|6l[| . To the 
best of our knowledge, however, in previous work viola- 
tions of dynamical scaling have been seen only in some 
systems with vector and more complex order parameter. 
In contrast, Lotka-Volterra models can be interpreted as 
systems with scalar order parameter, although the num- 
ber of equilibrium states N generally exceeds two, the 
characteristic value for Ising-type systems. 

Finally we note that presence of only two length scales 
exemplifies the mildest violation of classical single-size 
scaling. Generally, if scaling is violated one expects the 
appearance of an infinite number of independent scales, 
i.e., multiscaling p4] , |62| . Surprisingly, we found no ev- 
idence of multiscaling. Similar two-length scaling has 
been observed in the simplest one-dimensional system 
with vector order parameter, namely in the XY model 
j58| , and in the single-species annihilation with combined 
diffusive and convective transport (3^] . Indications of the 
three-length dynamical scaling have been reported in the 
context of coarsening Q and chemical kinetics |33]]|]. 

VIII. SUMMARY 

In this study, we addressed the dynamics of compet- 
itive immobile species forming a cyclic food chain. Wc 
first examined a cyclic model with asymmetric rules and 
symmetric initial conditions and have observed a drastic 
difference between the two extremes, corresponding to 
the complete graph ("infinite-dimensional") and to one- 
dimensional substrates. In the latter case, spatial inho- 
mogeneities develop, and the resulting kinetic behavior 
is very sensitive to the number of species. For a suffi- 
ciently small number of species, the system coarsens and 
is described by a set of exponents summarized in Table 
1. These exponents depend on the number of species and 
on the type of dynamics (parallel or sequential). Thus, 
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to describe coarsening in systems with non- conservative 
dynamics it is necessary to specify the details of the dy- 
namics. 

The time distribution of the number of mutations has 
also been investigated and we presented scaling argu- 
ments as well as an exact result for a particular case. 
We also treated symmetric interaction rules. This sys- 
tem is especially interesting when N = 4 as it provides 
a clear realization of the recently introduced notion of 
"persistent" spins in terms of the stationary interfaces. 
Finally, we discussed a relationship to the alternative 
reaction-diffusion equations description. While for the 
cyclic version both the lattice and the reaction-diffusion 
approaches have been found to be closely related, for 
the symmetric version very different results have emerged 
and a relationship with extremal dynamics has been es- 
tablished. 
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